function [pv,pIv,vv,bfail] =SuffNSolBellstf(PI,M,S,vbarv,bet,gam) 
% Bellman function iteration N-state model with fx rate 

N = max(size(vbarv));
bfail=1;
pbarv=(eye(N)-bet*gam*PI)\vbarv;

PIM =((PI.*M)*gam.*S);
pI = PIM*pbarv;
pbarU = pbarv;
palg = max(pI,pbarU);

k=1000000;

critm1 = 10e-9;
%Palg=zeros(N,k);
for j=1:k
    pI = PIM*palg;
    pbarU = vbarv+ bet*gam*PI*palg;
    palgn = max(pI,pbarU);
    if j>500
        diffm1=max(abs(palgn-palg));
        if diffm1<critm1
            'Bellman equation successfully converged'
            j
            palg=palgn;
            bfail=0;
            break
        end
    end
    palg=palgn;
end

pv=palg;
pIv = pI;
vv = max(vbarv,(((PI.*M)*gam.*S)-(bet*gam)*PI)*pv);

end

